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EXECUTIVE  SUMMARY 


This  is  the  final  report  for  SERDP  SEED  Project  MM-1642  and  it  covers  the 
research  results  and  developments  for  enhancement  of  magnetic  data  by  stable  down¬ 
ward  continuation  for  UXO  applications.  The  magnetic  method  is  one  of  the  two 
most  effective  geophysical  techniques  currently  in  use  for  UXO  detection  and  discrim¬ 
ination.  These  data  are  always  acquired  with  sensors  at  some  height  above  ground 
surface,  and  it  is  well  known  that  height  significantly  affects  magnetic  anomaly  shape, 
amplitude,  and  spatial  extent.  Consequently,  anomalies  due  to  multiple  metallic  tar¬ 
gets  may  overlap  at  a  given  height  above  the  ground  surface,  and  the  acquisition  noise 
may  significantly  decrease  the  signal-to-noise  ratio  (SNR)  of  data.  These  adverse  ef¬ 
fects  ultimately  mask  the  true  level  of  contamination  at  a  site  during  initial  wide  area 
assessment  (WAA),  as  well  decrease  overall  effectiveness  of  discrimination  during  the 
active  clearance  stage.  The  overall  objective  of  the  current  project  has  therefore  been 
to  develop  a  robust  algorithm  for  stable  downward  continuation  (SDC)  of  magnetic 
data  acquired  at  some  height  above  the  ground  to  reconstruct  the  magnetic  data  with 
a  higher  resolution  at  the  ground  surface. 

The  research  performed  over  the  first  year  of  the  project  progressed  satisfactorily 
with  the  major  research  tasks  completed  as  scheduled  and  a  working  SDC  algorithm 
developed  and  tested  on  synthetic  and  field  data.  The  algorithm  formulates  the 
downward  continuation  as  an  inverse  problem  using  Tikhonov  regularization  and  has 
the  flexibility  of  incorporating  the  expected  power  spectrum  of  UXO  anomalies.  The 
degree  of  regularization  is  estimated  automatically  using  the  well-established  methods 
in  linear  inverse  problems.  Applications  show  that  the  algorithm  can  reliably  estimate 
the  noise  in  the  data  and  reconstruct  the  magnetic  anomaly  at  ground  surface  within 
the  limitation  imposed  by  the  noise.  The  reconstructed  field  at  the  ground  surface 
exhibits  significant  enhancement  compared  to  the  original  data. 

iii 


The  second  year  of  the  project  significantly  advanced  upon  the  original  proposed 
research  by  incorporating  an  initial  set  of  inversion  algorithms  to  estimate  the  req¬ 
uisite  ensemble  depths  within  the  data’s  radial  power  spectrum  for  improved  SDC. 
In  particular,  development  of  recursive  Quenched  Simulated  Annealing  (QSA)  from 
SERDP  Project  MM-1638  has  been  adapted  for  this  purpose.  The  QSA  algorithm 
has  demonstrated  a  reliable  tool  for  modeling  the  radial  power  spectrum  of  UXO 
magnetic  data  using  non-linear  parametric  inversion.  As  a  result,  we  can  now  auto¬ 
matically  define  a  model  objective  function  for  SDC  to  quantify  the  conformity  of  the 
enhanced  data  with  expected  spectral  properties. 

In  this  report,  we  present  the  theoretical  and  practical  aspects  for  development  of 
stable  downward  continuation  in  UXO  applications.  We  provide  a  means  of  estimating 
one’s  expected  power  spectrum  within  their  UXO  data  through  parametric  inversion, 
and  demonstrate  successful  application  of  SDC  to  synthetic  and  field  data. 
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CHAPTER  1 


INTRODUCTION 


1.1  SERDP  Relevance 


This  project  addresses  the  need  for  signal  processing  procedures  as  outlined  in 
the  Statement  of  Needs  MMSEED-08-01,  and  seeks  to  develop  a  new  technology  for 
enhancing  magnetic  data  by  increasing  data  resolution  and  by  estimating  data  noise 
for  use  in  subsequent  analyses. 

The  magnetic  method  is  one  of  the  two  most  effective  geophysical  techniques 
currently  in  use  for  UXO  detection  and  discrimination.  The  use  of  magnetic  data 
in  UXO  applications  has  been  well  studied  and  there  are  well-documented  successful 
examples  of  locating  and  discriminating  UXO.  Detection  of  UXO  through  magnetic 
data  depends  on  identifying  dipole-like  anomalies,  while  discrimination  depends  upon 
inversion  to  recover  reliably  the  magnetic  source  parameters.  These  magnetic  data 
are  acquired  with  sensors  at  some  height  above  ground  surface.  It  is  well  known  that 
height  significantly  affects  magnetic  anomaly  shape,  amplitude,  and  spatial  extent. 
Consequently,  anomalies  due  to  multiple  metallic  targets  may  overlap  at  a  given 
height  above  the  ground  surface,  and  the  acquisition  noise  may  significantly  decrease 
the  signal-to-noise  ratio  (SNR)  of  data.  An  example  from  Kirtland  Air  Force  Base, 
NM  (Billings  et  ah,  2009),  is  presented  in  Figure  1.1.  The  data  in  this  example  are 
airborne  magnetic  data  collected  approximately  1  meter  above  the  ground  surface  over 
a  281  m  by  281  m  grid  and  it  demonstrates  such  overlapping  anomalies  and  decreased 
SNR.  These  adverse  effects  will  decrease  the  effectiveness  of  current  detection  and 
discrimination  methods  and  may  hamper  the  cost  reduction  of  UXO  clearance. 
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Figure  1.1.  Airborne  magnetic  data  collected  from  UXO  site  at  Kirtland  Air  Force 
Base,  NM.  The  data  are  in  nano  Tesla  and  it  illustrates  numerous  overlapping  anoma¬ 
lies  from  potential  UXOs.  The  data  clearly  demonstrates  the  need  for  enhancing 
airborne  magnetic  data  to  suppress  acquisition  noise  and  separate  overlapping  anoma¬ 
lies. 


1.2  Technical  Objective 

The  purpose  of  this  study  is  to  develop  a  robust  algorithm  for  stable  downward 
continuation  of  magnetic  data  acquired  at  some  height  above  the  ground  to  recon¬ 
struct  the  magnetic  data  with  a  higher  resolution  at  the  ground  surface,  and  also  to 
characterize  the  noise  in  ground-based  and  airborne  data.  The  specific  objectives  are: 

1.  Developing  a  robust  algorithm  for  stable  downward  continuation  of  magnetic 
data; 

2.  Feasibility  studies  of  application  to  airborne  magnetic  data,  especially  in  wide 
area  assessment  (WAA)  applications; 

3.  Feasibility  studies  of  enhancing  ground-based  magnetic  data;  and 

4.  Noise  estimation  and  its  use  in  improved  inversion  algorithms. 
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1.3  Project  Summary 


In  this  report,  we  summarize  the  major  research  accomplishments  of  SERDP 
MM-1642  SEED  Project.  In  short,  we  have  accomplished  the  main  objectives  of  this 
project,  and  we  have  made  significant  additional  contributions  to  advance  beyond  the 
initial  project  goals. 

In  Chapter  2,  we  present  the  main  development  of  stable  downward  continuation 
(SDC)  for  enhancing  magnetic  data  in  UXO  applications.  The  algorithm  is  formulated 
as  an  inverse  problem  using  a  regularization  approach,  and  it  has  been  developed  to 
incorporate  the  expected  power  spectrum  of  UXO  anomalies.  We  expand  beyond  the 
original  research  goals  of  MM-1642  in  Chapters  3.  We  have  identified  that  there  are 
distinct  ensembles  in  all  UXO  magnetic  data,  and  that  we  can  efficiently  recover  the 
predicted  power  spectrum  automatically  through  non-linear  parametric  inversion  for 
these  ensemble  based  models.  This  understanding  has  both  theoretical  implications 
and  practical  utility  in  processing  UXO  magnetic  data  and  extracting  anomalies. 
We  have  incorporated  latest  developments  in  Quenched  Simulated  Annealing  from 
SERDP  Project  MM-1638  for  this  stage  of  of  inverse  problem. 

In  Chapter  4,  we  identify  the  parameters  under  which  SDC  may  be  able  to 
increase  the  data  resolution  and  gain  an  understanding  of  the  limitations  of  the  algo¬ 
rithm.  We  implement  the  algorithm  on  data  sets  to  examine  the  data  enhancement 
results,  and  to  illustrate  the  possibility  of  using  the  algorithm  to  assess  the  noise  in 
the  data.  Chapter  5  closes  the  research  project  with  a  series  of  large-scale  field  exam¬ 
ples,  demonstrating  the  effectiveness  of  SDC  for  enhancing  airborne  magnetic  data 
from  the  ESTCP  wide  area  assessment  (WAA)  demonstration  project  (MM-0741)  at 
Kirtland  Air  Force  Base. 

We  conclude  the  report  in  Chapters  6  and  7  with  a  discussion  of  the  project 
developments,  results,  recommendations  and  publications/presentations  generated 
throughout  the  lifespan  of  SERDP  SEED  Project  MM-1642. 
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CHAPTER  2 


STABLE  DOWNWARD  CONTINUATION  (SDC) 


2.1  Introduction 

In  this  chapter,  we  present  an  algorithm  for  enhancing  magnetic  data  in  UXO 
applications  using  a  stable  downward  continuation  (SDC)  method.  The  algorithm 
formulates  the  downward  continuation  as  an  inverse  problem  using  Tikhonov  regu¬ 
larization  (Tikhonov  and  Arsenin,  1977)  and  has  the  flexibility  of  incorporating  the 
expected  power  spectrum  of  UXO  anomalies.  The  degree  of  regularization  is  esti¬ 
mated  automatically  using  the  well-established  methods  in  linear  inverse  problems. 
Numerical  tests  show  that  the  algorithm  can  reliably  estimate  the  noise  in  the  data 
and  reconstruct  the  magnetic  anomaly  at  ground  surface  within  the  limitation  im¬ 
posed  by  the  noise.  The  reconstructed  field  at  the  ground  surface  exhibits  significant 
enhancement  compared  to  the  original  data.  This  chapter  discusses  the  stable  down¬ 
ward  continuation  algorithm,  the  estimation  of  expected  power  spectrum,  the  choice 
of  regularization  level,  and  demonstrates  the  method  with  a  synthetic  data  set. 


2.2  SDC  Algorithm 


The  magnetic  data  at  two  observation  heights  are  related  by  the  upward  contin¬ 
uation  operation  (Blakely,  1996), 


Th(x,  y,  A  ft) 


!_  r  r  mxW)Ah  , 

2>r  J-o o  J- co  \(x'  -  xf  +  (</  -  yf  +  AVp/2  X  V 


(2.1) 


where  T0(x,y )  and  Th(x,y,h)  are  respectively  the  magnetic  data  at  two  observation 
heights  separated  by  a  vertical  distance  Ah.  Applying  a  two-dimensional  Fourier 
transform  to  Equation  (2.1)  yields  a  simpler  form  in  which  the  Fourier  transforms  of 
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the  two  quantities  are  related  to  each  other  by  a  simple  upward  continuation  operator 


Th(ux,Uy,Ah)  =  e  AhuJrf0(ux,uy),  (2.2) 

where  T0(u>x,ujy )  denotes  the  Fourier  transform  of  T0(x,y )  ,  ujx  and  uy  are  wavenum¬ 
bers  in  x—  and  y— direction  respectively,  and  ixr  =  ^JoJx  +  0Jy  is  the  radial  wavenum¬ 
ber.  The  upward  continuation  operator  attenuates  the  high-frequency  content  of  a 
magnetic  anomaly  as  a  function  of  height.  These  effects  can  be  partially  reversed 
by  performing  downward  continuation,  which  is  the  inverse  of  Equation  (2.2),  to  re¬ 
construct  the  field  at  the  ground  surface.  The  ground  surface  is  the  lowest  possible 
observation  height. 

Although  the  operation  of  downward  continuation  is  mathematically  valid  in  the 
source  free  region  above  the  ground  surface  (Blakely,  1996),  it  is  numerically  unstable 
because  of  the  presence  of  high-frequency  noise,  ffowever,  we  can  stabilize  the  process 
by  formulating  it  as  an  inverse  problem  and  derive  a  regularized  operator  (Huestis  and 
Parker,  1979).  We  will  carry  out  this  operation  in  the  wavenumber  domain  because 
of  the  numerical  efficiency. 

Let  us  treat  the  Fourier  transform  Th(u>x,ujy,  Ah)  of  observed  magnetic  anomaly 
at  height  Ah  as  the  data,  and  the  sought  Fourier  transform  of  anomaly  at  the  ground 
surface  as  the  model.  The  upward  continuation  operator  e~Ahuir  then  defines  the  for¬ 
ward  mapping.  The  inverse  problem  is  stated  as  hireling  a  well-behaved  T0(ujx,  c oy)  that 
has  certain  spectral  characteristics  and  reasonably  reproduces  the  data  Th(u;x,  coy ,  Ah) 
by  Equation  (2.2).  Assuming  the  magnetic  data  that  would  be  observed  at  the  ground 
surface  has  an  expected  power  spectrum  of  Po(cox,  ojv)  ,  we  can  define  a  model  objec¬ 
tive  function  (j)rn  to  quantify  the  conformity  of  the  model  with  the  expected  spectral 
property, 

(pm  j Pq (^-T ,  )  | Tq (ujx ,  coy )  |  dioxduy ,  (2.3) 

where  the  reciprocal  of  the  expected  power  spectrum  provides  the  weighting  function. 


6 


An  assumed  P0(ujx,ujy )  that  decays  with  increasing  wavenumber  corresponds  to  an 
increased  weight  at  higher  wavenumber  for  the  final  model.  Consequently,  the  final 
model  from  the  inversion  will  have  the  similar  spectral  decay.  Such  a  spectral  decay 
corresponds  to  a  smooth  magnetic  held  in  the  space  domain.  Thus,  the  use  of  inverse 
spectral  decay  as  a  weighting  function  provides  a  versatile  approach  to  incorporate 
prior  information  about  the  magnetic  data.  We  will  discuss  this  aspect  in  more  detail 
in  the  next  section. 

We  define  a  data  misfit  function  as, 

4>d=  I  I \T0(ujx,ujy)e-AhuJr  -  T0(ujx,ujy,h)\2  duxduy,  (2.4) 

which  is  in  effect  the  scaled  version  of  data  misfit  in  the  space  domain  by  Parseval’s 
theorem.  The  lack  of  a  normalization  factor  in  this  definition  also  implies  the  as¬ 
sumption  of  independent  noise  in  the  space  domain,  which  translates  to  a  constant 
noise  level  across  the  wavenumber  band  of  the  data. 

Utilizing  the  formalism  of  Tikhonov  regularization,  the  inverse  problem  of  recon¬ 
structing  the  magnetic  anomaly  at  the  ground  surface  then  becomes  one  of  optimiza¬ 
tion: 

min  -.(j)  =  (j)d  +  //0m,  (2.5) 

where  /i  is  the  regularization  parameter  and  it  determines  the  trade-off  between  fitting 
the  data  and  the  conformity  of  model  with  the  expected  spectral  properties.  Solving 
the  minimization  with  the  optimal  value  of  regularization  parameter  /j  produces  the 
desired  Fourier  transform  of  downward  continued  data. 

The  solution  of  Equation  (2.5)  involves  both  complex  data  and  a  complex  model. 
Therefore,  we  resort  to  the  calculus  of  variation  to  carry  out  the  formal  minimization. 
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The  result  is  a  simple  decoupled  system  and  we  have  a  closed-form  inverse  solution: 


DAhuiR 


1  +  /iP0(cox,uJy)e2AhujR 


Th  i  ^ y )  AH) . 


(2.6) 


Performing  a  2D  inverse  Fourier  transform  of  To(ujx,cjy)  then  yields  the  sought  mag¬ 
netic  data  at  the  ground  surface. 


The  choice  of  regularization  parameter  /i  determines  the  level  of  data  fit  and 
therefore  the  smoothness  of  the  downward  continued  field.  If  we  have  precise  knowl¬ 
edge  of  the  statistics  of  errors  in  the  data,  we  can  define  the  expected  value  of  the 
data  misfit  in  Equation  (2.4)  and  the  optimal  regularization  parameter  must  be  the 
one  that  yields  the  expected  misfit.  In  general,  we  have  little  precise  information 
about  the  noise  characteristics,  and  this  approach  is  not  applicable.  In  such  cases, 
we  must  resort  to  other  methods  to  determine  /i.  The  process  then  is  equivalent  to 
estimating  the  noise  level  in  the  data. 


There  are  several  approaches  for  performing  the  estimation.  We  apply  the  L- 
curve  criterion  (Hansen,  1992).  The  L-curve  criterion  is  a  heuristic  method  based 
on  the  behavior  of  Tikhonov  curve.  It  was  observed  that,  when  plotted  on  a  log-log 
scale,  the  data  misfit  as  a  function  of  model  objective  function  exhibits  a  character¬ 
istic  corner.  As  the  degree  of  regularization  decreases  towards  this  corner  point,  the 
model  objective  function  changes  very  little  while  the  misfit  is  being  reduced  greatly. 
Further  decrease  in  the  degree  of  regularization  beyond  this  point  would  result  in 
rapid  increase  in  the  model  complexity  with  little  reduction  in  the  data  misfit.  The 
L-curve  criterion  considers  the  corner  point  to  be  the  best  compromise  between  ex¬ 
tracting  maximum  amount  of  signal  and  being  least  affected  by  noise.  It  therefore 
states  that  the  optimal  solution  for  a  given  regularized  inverse  problem  is  the  one 
that  corresponds  to  the  corner  point  of  the  Tikhonov  curve.  The  corner  is  defined 
numerically  as  the  point  with  the  maximum  curvature. 


2.3  Power  spectrum  of  magnetic  data  in  UXO 


The  formulation  presented  in  the  preceding  section  is  general,  and  one  can  poten¬ 
tially  use  any  assumed  form  of  power  spectrum  to  define  the  model  objective  function 
in  order  to  stabilize  the  inversion.  The  weighting  function  based  the  reciprocal  of  as¬ 
sumed  power  spectrum  has  the  ultimate  effect  of  re-shaping  the  spectral  content  of 
the  downward  continued  magnetic  data.  The  faster  decaying  power  spectrum  corre¬ 
sponds  to  a  smoother  final  result.  For  example,  assuming  P0  =  js  equivalent 

to  requiring  the  downward  continued  data  to  be  smooth  by  the  first-order  derivative 
in  the  space  domain.  In  the  absence  of  specific  information  regarding  the  nature  of 
the  source,  such  a  generic  assumption  is  useful  for  obtaining  first  order  results.  To 
obtain  the  best  possible  result,  however,  one  must  incorporate  the  correct  spectral 
property  consistent  with  the  data  set  in  question. 

For  our  current  application,  we  choose  to  use  the  theoretical  radially  averaged 
power  spectrum.  Consider  a  set  of  UXO  items  producing  the  magnetic  data  in  ques¬ 
tion.  Each  UXO  is  represented  by  a  magnetic  dipole  to  a  high  degree  of  accuracy. 
Thus,  we  have  an  ensemble  of  random  dipoles  with  varying  dipole  moments  and  burial 
depths.  Spector  and  Grant  (1970)  showed  elegantly  that,  within  a  scaling  factor,  the 
radially  averaged  power  spectrum  of  the  magnetic  data  produced  by  an  ensemble  is 
the  same  as  that  of  the  data  produced  by  an  ’’average”  member  of  the  ensemble.  For 
a  dipole  source,  this  power  spectrum  has  the  following  form, 

P(ur)  oc  Aco2e~2huJ\  (2.7) 

where  h  denotes  the  vertical  distance  between  the  ensemble  average  and  the  obser¬ 
vation  plane  and  A  is  the  scaling  factor  depending  on  the  source  strength.  It  follows 
naturally  that  we  use  Equation  (2.7)  to  define  the  spectral  weighting  for  enhancing 
magnetic  data  in  UXO  applications. 

For  practical  applications,  one  must  estimate  the  depth  of  the  ensemble  average 
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below  the  observation  plane.  This  can  be  accomplished  by  performing  a  least-squares 
fit  between  Equation  (2.7)  and  the  radially  averaged  power  spectrum  of  the  data. 
To  account  for  the  presence  of  noise,  we  may  also  include  a  constant  noise  power  to 
Equation  (2.7).  In  general,  the  magnetic  data  are  more  complex  than  those  produced 
by  a  single  ensemble.  For  example,  the  magnetic  soil  will  manifest  itself  in  the  data 
and  exhibit  as  an  ensemble  at  a  depth  near  the  ground  surface.  There  could  also  be 
deeper  geology  at  much  greater  depths  below  the  UXOs.  Thus,  the  power  spectrum  of 
a  data  set  is  typically  consistent  with  that  of  multiple  source  ensembles,  each  having 
a  different  depth  and  source  strengths.  In  addition,  there  is  usually  a  random  noise 
component.  To  obtain  the  correct  power  spectrum  for  use  in  the  stable  downward 
continuation  algorithm,  we  must  estimate  these  parameters  numerically. 

Once  the  average  depth  h  of  the  UXO  ensemble  is  obtained  through  this  process, 
the  power  spectrum  P0(uj’xujy )  used  in  defining  the  model  objective  function  is  given 
by 


P0(ux,uy)=u2re-2^-Ah^j  (2.8) 

where  the  vertical  distance  h  —  Ah  defines  the  ensemble  depth  below  the  ground 
surface  to  which  we  downward  continue  the  data. 

As  an  example,  Figure  2.1  illustrates  the  power  spectral  estimation  using  a  field 
data  set  from  Kirtland  Air  force  Base,  NM.  The  blue  line  is  the  actual  power  spectrum 
calculated  from  the  FFT  of  the  magnetic  data  over  a  30  m  by  30  m  area.  The  red 
line  is  the  radially  averaged  power  spectrum  consisting  of  two  distinct  ensembles  and 
a  constant  noise  power.  The  first  ensemble’s  average  depth  is  consistent  with  the 
buried  UXO  in  the  area,  and  the  second  ensemble  corresponds  to  the  observation 
height  above  the  ground.  The  agreement  between  the  actual  and  theoretical  power 
spectra  with  two  ensembles  is  remarkable. 
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Figure  2.1.  Radially  averaged  power  spectrum  and  its  parametric  representation.  The 
blue  curve  is  the  numerically  calculated  power  spectrum  of  a  subset  of  magnetic  data 
at  Kirtland.  The  red  line  represents  the  radially  averaged  power  spectrum  and  it 
consists  of  three  distinct  ensemble  averages  plus  a  constant  noise  power 

2.4  Synthetic  example 

We  now  proceed  to  illustrate  the  algorithm  using  a  synthetic  example.  We  first 
illustrate  the  details  of  the  algorithm  and  show  the  effectiveness  of  the  stable  down¬ 
ward  continuation.  We  then  use  this  example  to  evaluate  the  range  of  observation 
height  from  which  the  enhancement  is  viable. 

We  use  a  model  consisting  of  two  buried  dipoles  separated  horizontally  by  2  m 
and  buried  at  the  same  depth  of  0.5  m  below  the  surface.  Assuming  an  observation 
height  of  2  m  above  the  ground  surface  and  adding  0.5  nT  of  Gaussian  noise,  we 
simulated  the  observed  data  with  0.1-m  spacing  shown  in  Figure  2.2.  At  this  height, 
the  magnetic  data  shows  one  broad  anomaly.  In  UXO  applications,  this  might  indicate 
the  presence  of  clutters,  but  it  would  be  difficult  to  discern  how  many  buried  targets 
are  present. 

Estimating  a  radial  power  spectrum  and  then  applying  the  stable  downward  con¬ 
tinuation  yields  the  result  shown  in  the  same  figure.  The  continued  data  are  smooth 
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Figure  2.2.  Results  of  stable  downward  continuation  applied  to  magnetic  data  gen¬ 
erated  for  two  closely  spaced  dipoles.  The  upper  left  panel  is  the  original  data  with 
noise  added.  The  upper  right  panel  is  the  enhanced  data  from  SDC  processing.  The 
lower  left  image  is  the  predicted  data  at  the  original  observation  height,  and  it  is  a 
good  fit  to  the  original  noisy  data.  The  lower  right  panel  is  the  noise  estimate  by 
subtracting  the  predicted  and  original  data. 


Figure  2.3.  The  top  panels  show  respectively  the  data  misfit  ((f>d)  and  model  objective 
function  ((f) m)  for  the  downward  continuation  shown  in  Figure  2.2.  The  Tikhonov 
curve  in  the  lower-left  panel  has  a  well-defined  corner  that  is  located  by  the  peak 
curvature  point. 
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Figure  2.4.  Comparison  between  the  true  anomaly  at  the  ground  surface  and  the 
stably  downward  continued  anomaly  for  the  example  shown  in  Figure  2.2. 

and  well  behaved.  Two  expected  anomalies  are  clearly  visible.  For  comparison,  we 
also  show  the  predicted  data,  which  is  a  de-noised  version  of  the  observed  data.  The 
estimated  noise  given  by  the  difference  between  the  observed  and  predicted  data  is 
representative  of  the  statistics  of  the  actual  noise. 

Figure  2.3  shows  more  details  regarding  the  performance  of  the  algorithm.  To 
search  for  an  optimal  regularization  parameter,  we  have  tested  a  range  of  values.  The 
data  misfit  exhibits  the  typical  monotonic  increase  with  the  regularization  parame¬ 
ter,  while  the  corresponding  model  objective  function  decreases  monotonically.  The 
Tikhonov  curve  has  a  well  defined  corner,  as  indicated  by  a  single  peak  in  curvature. 
This  optimal  regularization  level  yielded  the  result  in  Figure  2.2. 

As  a  final  assessment  in  this  example,  we  compare  the  true  anomaly  from  the  two 
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dipoles  at  the  ground  surface  with  the  downward  continued  data  in  Figure  2.4.  We 
note  that  the  downward  continued  data  are  smoother  and  the  amplitude  is  greatly 
reduced.  This  is  expected  since  the  continuation  process  must  attenuate  a  great  deal 
of  high-frequency  content  to  obtain  a  stable  solution.  Despite  the  smooth  appearance 
and  low  amplitude,  the  enhancement  is  significant  since  the  two  anomalies  are  clearly 
visible. 

2.5  Summary 

We  have  presented  a  stable  downward  continuation  method  for  enhancing  mag¬ 
netic  data  acquired  for  UXO  detection  and  discrimination.  The  key  components  are 
the  expected  power  spectrum  and  automatic  estimation  of  noise  level  in  the  data. 
The  algorithm  formulates  the  downward  continuation  as  an  inverse  problem  using 
Tikhonov  regularization  and  has  the  flexibility  of  incorporating  the  expected  power 
spectrum  of  UXO  anomalies.  The  degree  of  regularization  is  estimated  automatically 
using  the  well-established  methods  in  linear  inverse  problems.  Numerical  tests  show 
that  the  algorithm  can  reliably  estimate  the  noise  in  the  data  and  reconstruct  the 
magnetic  anomaly  at  ground  surface  within  the  limitation  imposed  by  the  noise. 

In  the  next  chapter,  we  present  advances  to  the  project  beyond  the  proposed 
research,  by  identifying  an  ensemble  based  representation  for  UXO  magnetic  data,  and 
by  modifying  recent  inversion  technologies  to  identify  the  expected  spectral  properties 
of  this  model  representation  as  necessary  for  SDC. 
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CHAPTER  3 


POWER  SPECTRUM  ESTIMATION 


In  the  previous  chapter,  we  demonstrated  the  approach  of  stable  downward  con¬ 
tinuation  to  enhance  magnetic  data  in  UXO  applications.  The  basic  algorithm  was 
developed  and  tested  throughout  the  first  year  of  the  project.  The  two  key  compo¬ 
nents  are  the  expected  power  spectrum  and  automatic  estimation  of  noise  level  in 
the  data.  The  second  year  of  the  project  significantly  advanced  upon  the  original 
proposed  research  by  incorporating  an  initial  set  of  inversion  algorithms  to  estimate 
the  requisite  ensemble  depths  and  noise  level  within  the  data’s  radial  power  spectrum 
for  improved  SDC.  In  particular,  development  of  recursive  Quenched  Simulated  An¬ 
nealing  (QSA)  from  SERDP  Project  MM-1638  (Krahenbuhl  et  ah,  2010)  has  been 
adapted  for  this  purpose.  The  QSA  algorithm  has  demonstrated  a  reliable  tool  for 
modeling  the  radial  power  spectrum  of  UXO  magnetic  data  using  non-linear  para¬ 
metric  inversion.  As  a  result,  we  can  now  automatically  define  a  model  objective 
function  for  SDC  to  quantify  the  conformity  of  the  enhanced  data  with  expected 
spectral  properties.  The  procedure  can  be  readily  transferred  for  use  in  practice. 

In  this  chapter,  we  provide  a  brief  background  on  the  power  spectrum  estimation 
problem,  its  formulation  as  a  parametric  inverse  problem,  solution  by  QSA.  We  also 
discuss  the  underlying  ensemble  models  we  have  adopted  in  representing  the  power 
spectra  of  UXO  magnetic  data  and  its  utility  in  understanding  the  source  components 
of  magnetic  data.  Although  the  estimation  of  the  power  spectrum  is  a  basic  com¬ 
ponent  of  the  proposed  work,  what  we  have  accomplished  is  a  significant  addition. 
The  understanding  achieved  through  the  ensemble  representation  of  magnetic  data 
has  both  theoretical  implications  and  practical  utility  in  processing  UXO  magnetic 
data  and  extracting  anomalies.  We  have  observed  that,  in  all  UXO  magnetic  data  we 
have  examined,  there  are  invariably  distinct  ensembles  of  magnetic  sources  present 
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and  each  ensemble  has  a  good  correspondence  with  a  source  group  such  as  magnetic 
soil,  UXOs,  and  underlying  geology.  This  observation  and  understanding  provide  a 
practical  and  informative  approach  to  examining  magnetic  data  and  separating  signal 
components  in  UXO  applications. 


3.1  Introduction 

For  our  current  inverse  problem,  we  have  opted  to  use  the  theoretical  radially 
averaged  power  spectrum.  Consider  a  set  of  UXO  items  producing  the  magnetic 
data  in  question.  Each  UXO  is  represented  by  a  magnetic  dipole  to  a  high  degree  of 
accuracy.  Thus,  we  have  an  ensemble  of  random  dipoles  with  varying  dipole  moments 
and  burial  depths.  Spector  and  Grant  (1970)  showed  elegantly  that,  within  a  scaling 
factor,  the  radially  averaged  power  spectrum  of  the  magnetic  data  produced  by  an 
ensemble  is  the  same  as  that  of  the  data  produced  by  an  ”  average”  member  of  the 
ensemble.  For  a  dipole  source,  this  power  spectrum  has  the  following  form, 

P{ur)  oc  (Auj2re-2h“r)  +  Pn7  (3.1) 

where  h  denotes  the  vertical  distance  between  the  ensemble  average  and  the  observa¬ 
tion  plane,  A  is  the  scaling  factor  depending  on  the  source  strength,  ujr  is  the  radial 
wavenumber,  and  P/v  is  a  constant  noise  power.  It  follows  naturally  that  we  use 
Equation  (3.1)  to  define  the  spectral  weighting  for  enhancing  magnetic  data  in  UXO 
applications. 

In  general,  one  would  expect  there  are  other  magnetic  sources  present  at  a  given 
UXO  site.  Correspondingly,  it  would  be  reasonable  to  use  multiple  ensembles  to 
represent  the  power  spectrum  of  a  data  set.  Indeed,  we  have  observed  this  to  be  true  in 
all  data  sets  we  have  examined  within  this  project  and  other  related  projects,  such  as 
MM-1638.  At  least  three  ensembles  are  present.  The  first  tends  to  be  associated  with 
the  soil  magnetism  and  metallic  objects  scatter  on  the  surface;  the  second  associated 
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with  UXO  and  other  metallic  objects  buried  in  the  ground;  and  the  third  with  geologic 
responses  at  larger  depth.  The  magnetic  sources  in  the  first  two  ensembles  are  depth 
limited  and  therefore  has  the  dipole  like  behavior.  Thus  the  form  in  Equation  (3.1) 
is  appropriate.  The  third  ensemble,  however,  appears  to  be  depth  unlimited  for  the 
scale  length  relevant  to  UXO  work.  The  corresponding  radial  power  spectrum  has 
the  form  of 

P{ur)  oc  Ae~2huJ\  (3.2) 

where  the  lack  of  the  factor  u2  describes  the  depth  unlimited  nature  of  the  sources. 

We  have  concluded  in  the  project  that  the  basic  radial  power  spectral  form  should 
consist  of  four  terms: 


P(ur)  =  A1u}?e~2hlUr  +  A2uj2e-2h2UJr  +  A3e~2h +  PN.  (3.3) 


Given  this  form,  we  structure  the  problem  as  a  non-linear  parametric  inverse  problem 
that  seeks  to  recover  the  seven  parameters  (A\,  A2,  A3,  hi,  h2,  h3,  and  Pn )  for  the  three 
ensemble  system.  The  associated  objective  function  for  the  minimization  problem  is 
then  expressed  as  the  2-norm  difference  between  the  logarithm  of  observed  ( dobs )  and 
predicted  radial  power  spectra(dpre): 


d: 


obs 


dpr 


(3.4) 


The  procedure  begins  by  performing  fast  Fourier  transform  (FFT)  of  the  gridded 
magnetic  data,  compute  the  periodogram  (square  of  the  modulus  of  the  FFT  )  as 
an  approximation  to  the  2D  power  spectrum,  and  then  then  compute  the  radially 
average  power  spectrum.  The  use  of  logarithm  is  necessary  since  radial  power  spectra 
varies  over  many  orders  of  magnitude. 

Given  the  exponential  and  logarithmic  functions  involved  in  the  objective  func¬ 
tion  in  Equation  (3.3),  this  is  a  highly  non-linear  parameter  estimation  problem. 
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The  solution  to  this  parametric  inverse  problem  is  nontrivial.  Although  the  stan¬ 
dard  Gauss-Newton  approach  (Nocedal  and  Wright,  1999)  works  to  some  extent,  it 
usually  requires  significant  user  intervention  in  the  form  of  initial  guess  to  ensure 
that  the  minimization  procedure  starts  within  the  convergence  radius  of  the  Gauss- 
Newton  approach.  To  overcome  this  difficulty,  we  modify  an  optimization  algorithm 
developed  through  SERDP  Project  MM-1638  for  scanning  the  solution  landscape. 
The  method  is  a  recursive  Quenched  Simulated  Annealing  (QSA)  (Krahenbuhl  et  ah, 
2010),  and  we  next  present  background  information  on  QSA  prior  to  demonstrating 
its  application  to  the  problem  at  hand. 

3.2  Quenched  Simulated  Annealing 

Simulated  Annealing  (SA)  is  global  search  technique  (Kirkpatrick,  1984)  well 
suited  for  non-linear  parametric  inversion  such  as  we  address  in  this  project.  SA 
formulation  is  designed  to  mimic  the  process  of  chemical  annealing,  where  the  final 
energy  state  of  a  crystal  lattice  is  determined  by  its  rate  of  cooling  through  the 
melting  point.  To  achieve  a  lower  energy  state  with  highly  ordered  crystals,  the 
material  is  usually  cooled  slowly.  When  the  material  is  cooled  too  rapidly,  the  lattice 
may  not  reach  the  lowest  possible  energy  state.  The  analogy  in  an  inversion  for  this 
latter  case  is  pre-mature  convergence  where  the  final  solution  has  missed  the  desired 
global  minimum.  For  geophysical  inversion,  the  SA  typically  starts  with  a  model  at 
random,  and  calculates  the  model’s  energy  based  on  its  objective  value  as  defined 
in  the  previous  chapter.  Perturbations  are  then  applied  to  the  model  and  the  new 
objective  values  are  calculated  at  each  iteration.  If  the  new  objective  value  decreases 
or  remains  the  same,  the  model  is  accepted  as  a  replacement.  If  the  objective  value 
increases,  the  model  is  accepted  by  a  thermally  controlled  probability  function  often 
referred  to  as  the  Metropolis  criterion: 


p  =  e(-A  */T)? 


(3.5) 
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where  A0  is  the  difference  between  objective  values  of  the  old  and  new  models, 
and  T  is  a  temperature  parameter  designed  to  decay  (or  cool)  over  time.  For  more 
information  on  SA,  the  reader  is  referred  to  Metropolis  et  ah  (1953);  Kirkpatrick 
(1984)  and  Nnlton  and  Salamon  (1988).  Applications  of  SA  specific  to  geophysical 
inversion  are  available  in  Sen  et  al.  (1995);  Nagihara  and  Hall  (2001);  Scales  et  al. 
(1992);  Sambridge  and  Mosegaard  (2002),  and  Martin  et  al.  (2007). 

We  have  opted  to  use  a  modified  SA,  called  Quenched  Simulated  Annealing 
(QSA),  as  a  local  search  tool  for  this  (Krahenbuhl  et  ah,  2010).  QSA  in  its  simplest 
form  is  Simulated  Annealing  described  previously,  without  the  Metropolis  or  any 
other  cooling  criteria  such  as  Equation  (3.5).  The  algorithm  works  as  follows.  After  a 
change  is  performed  to  a  model  by  QSA,  i.e.  one  parameter  is  perterbed,  the  objective 
value  of  the  new  solution  is  calculated  through  Equation  (3.4).  If  the  change  does 
not  increase  the  total  objective  value  the  perturbation  is  accepted.  Therefore,  the 
algorithm  only  accepts  downhill  and  lateral  moves.  Implemented  recursively  with 
random  starting  models,  the  method  has  a  high  likelihood  of  identifying  any  potential 
local  minima  within  the  solution  space,  as  well  as  the  desired  global  minimum. 

As  an  illustration,  we  present  the  parameter  fit  to  the  power  spectrum  of  a  held 
data  set.  Figure  3.1  shows  the  manual  fit  to  the  data  power  spectrum  (red  dots) 
using  the  manually  adjusted  parameters  (solid  line).  Note  that  we  can  reproduce  the 
general  shape  of  the  power  spectrum  through  this  manual  approach.  However,  as  we 
illustrate  next,  recovery  of  the  ensemble  parameters  in  Equation  (3.3)  through  robust 
inversion  yields  significant  improvement  over  the  manual  approach. 

Figure  3.2  demonstrates  improved  recovery  of  the  three  ensemble  parameters 
for  Kirtland  over  manual  fitting.  This  is  demonstrated  by  the  well-resolved  power 
spectrum  (solid-line)  in  comparison  to  the  true  data  (red-points).  The  high-quality 
fit  of  the  data  following  inversion  verifies  that  we  have  successfully  recovered  the 
expected  spectral  properties  for  the  Kirtland  AFB  data. 
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Observed  (r.);  Predicted  (b) 


Figure  3.1.  Manual  fitting  of  example  radial  power  spectrum  of  airborne  magnetic 
data.  The  red-points  are  the  true  data  calculated  for  the  site.  The  solid-line  is  the 
predicted  data  by  manual  picking  of  ensemble  parameters  prior  to  stable  downward 
continuation.  The  manual  approach  can  adequately  reproduce  the  general  trend  of 
the  data;  however,  a  better  fit  of  the  data  can  be  achieved  through  robust  inversion 
for  the  ensemble  parameters. 


Observed  (r.) ;  Predicted  (b) 


Figure  3.2.  QSA  fitting  of  example  radial  power  spectrum  of  airborne  magnetic  data. 
The  red-points  are  the  true  data  calculated  for  the  site.  The  solid-line  is  the  predicted 
data  generated  from  QSA  inversion  for  the  seven  unknown  parameters  in  equation 
(3.3).  Parametric  inversion  has  significantly  improved  recovery  of  the  three  ensemble 
system  over  the  manual  approach  from  Figure  3.1 
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3.3  Field  example 


We  now  present  the  result  of  power  spectrum  estimation  using  the  three  ensemble 
model  applied  to  airborne  magnetic  data  collected  for  UXO  remediation  efforts  at 
Kirtland  Air  Force  Base,  NM  (Billings  et  ah,  2009).  The  original  magnetic  data, 
collected  approximately  one  meter  above  the  ground  is  presented  in  Figure  3.3.  This 
is  a  subset  of  the  entire  data  set.  The  data  are  gridded  at  an  interval  0.25  m  given 
the  nominal  1-m  line  spacing. 

The  corresponding  radial  power  spectral  fitting  is  show  in  Figure  3.4.  There 
are  three  distinct  ensembles  seen  in  the  spectrum  and  a  rather  constant  noise  power 
is  clearly  seen  in  high-frequency  band.  We  note  that  the  gridding  interval  of  0.25 
m  could  be  slightly  smaller  than  the  data  requires,  because  the  fluctuation  in  the 
spectrum  in  the  highest  frequency  band  is  clearly  numerical  noise. 

The  depths  of  the  three  ensembles  below  the  data  surface  are  respectively:  0.7, 
2.5,  and  11.5  m.  We  have  interpreted  these  results  as  follows.  The  ensemble  with  a 
depth  of  0.7  m  corresponds  to  general  noise  in  the  survey,  which  consists  of  surhcial 
variation  of  soil  magnetism  and  clutters.  Another  source  of  the  noise  could  be  the 
altitude  variation  and  positioning  of  the  airborne  platform.  An  average  depth  below 
the  observation  of  0.7  m,  which  is  slightly  above  the  ground,  is  quite  reasonable.  The 
second  ensemble  has  an  average  depth  of  2.5  m.  Given  the  flight  height  of  1.0  to 
1.25  m,  this  ensemble  is  about  1.25  m  below  the  ground  surface.  This  ensemble  is 
precisely  the  expression  of  UXO  and  buried  metallic  objects  detected  in  this  survey. 
The  third  ensemble  has  an  apparent  depth  of  11.5  m  and  it  is  reasonable  to  identify 
this  group  as  geologic  in  origin. 

3.4  Summary 

In  this  section,  we  presented  a  method  for  automatically  and  reliably  estimating 
the  spectral  properties  of  UXO  magnetic  data.  The  approach  implements  a  recursive 
Quenched  Simulated  Annealing  (QSA)  leveraged  from  SERDP  Project  MM-1638, 
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Figure  3.3.  Airborne  magnetic  data  collected  over  UXO  site  at  Kirtland  Air  Force 
Base,  NM.  The  data  show  high  concentration  of  potential  UXOs  at  the  western  end 
of  the  sub-grid,  with  numerous  overlapping  anomalies. 


Figure  3.4.  Radial  power  spectrum  data  calculated  for  the  Kirtland  AFB  airborne 
magnetic  data  presented  in  Figure  3.3.  The  red-points  are  the  true  data  calculated 
for  the  site.  The  solid-line  is  the  predicted  data  generated  from  QSA  inversion  for  the 
seven  unknown  parameters  in  equation  (3.3).  Parametric  inversion  has  significantly 
improved  recovery  of  the  three  ensemble  system  over  the  manual  approach  from  Figure 
3.1 
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and  solves  a  non-linear  parametric  inverse  problem  to  recover  the  seven  parameters 
(including  depths  and  noise  power)  of  a  three  ensemble  system.  As  a  result  of  these 
developments,  we  can  now  automatically  define  a  model  objective  function  for  sta¬ 
ble  downward  continuation  (SDC)  to  quantify  the  conformity  of  the  enhanced  data 
with  expected  spectral  properties.  In  the  next  chapter,  we  apply  the  SDC  approach 
to  two  data  sets  in  order  to  understand  the  limitations  of  the  algorithm,  to  demon¬ 
strate  increased  data  resolution  through  SDC,  and  to  evaluate  SDC  as  an  method  for 
estimating  noise  in  UXO  magnetic  data. 
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CHAPTER  4 


ANOMALY  SEPARATION  AND  DATA  ENHANCEMENT 


In  this  chapter,  we  apply  the  stable  downward  continuation  algorithm  developed 
in  the  preceding  two  data  sets.  The  goals  are  to  understand  the  performance  of  the 
algorithm  and  to  evaluate  the  benefit  that  can  be  derived  from  the  application  of  such 
an  algorithm.  The  first  data  set  is  the  synthetic  data  from  the  two-dipole  model  used 
in  Chapter  2  and  we  will  use  it  to  investigate  the  parameters  under  which  downward 
continuation  may  be  able  to  increase  the  data  resolution  and  gain  an  understanding 
of  the  limitations  of  the  algorithm.  The  second  data  set  is  a  small  subset  from  the 
Kirtland  WAA  survey,  and  we  use  it  to  examine  the  data  enhancement  and  illustrate 
the  possibility  of  using  the  algorithm  to  assess  the  noise  in  the  data. 

4.1  Anomaly  Separation  and  Limitation  of  SDC 

We  first  return  to  the  two-dipole  model  used  previously  and  investigate  the  limit 
of  the  observation  height  from  which  the  stable  downward  continuation  can  separate 
adjacent  dipole  anomalies.  The  ability  to  separate  and  distinguish  adjacent  anomalies 
ultimately  determines  the  increased  effective  resolution  of  magnetic  data  for  UXO 
detection. 

We  have  carried  out  the  simulation  for  a  number  of  different  observation  heights 
in  order  to  examine  this  aspect.  The  criterion  is  the  ability  to  identify  the  presence 
of  two  separate  anomalies  from  simulated  observations  at  different  heights.  The  two 
dipoles  are  buried  at  the  depth  of  0.5  m  below  the  ground  and  separated  by  2.0  m 
horizontally.  We  simulate  data  at  observation  heights  ranging  from  1  m  to  3  m  at  an 
interval  of  0.25  m.  Gaussian  random  noise  was  added  to  each  data  set.  The  data  sets 
are  then  stably  downward  continued  to  the  ground  surface  by  automatically  choosing 
the  optimal  regularization  parameter  in  each  case.  To  assess  the  quality  of  downward 
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Figure  4.1.  The  cross-correlation  coefficients  between  the  stably  downward  continued 
anomaly  and  the  true  anomaly  at  the  ground  surface. 


continued  anomaly,  we  calculate  its  cross-correlation  coefficient  between  the  true  mag¬ 
netic  anomaly  at  the  ground  surface.  We  also  visually  examine  the  continued  held 
for  the  presence  of  two  separate  anomalies.  Figure  4.1  displays  the  correlation  coef¬ 
ficient  as  a  function  of  the  observation  height.  Visual  inspection  indicates  that  the 
continued  anomaly  shows  two  separate  anomalies  at  a  height  equal  to,  or  less  than, 
2.5  m.  Correspondingly,  we  note  that  the  correlation  coefficient  continues  to  decrease 
for  observations  above  2.5  m.  This  set  of  simulations  therefore  shows  that  the  maxi¬ 
mum  observation  from  which  we  can  enhance  the  magnetic  data  by  stable  downward 
continuation  and  discern  two  separate  anomalies  is  2.5  m.  To  illustrate  this  limit, 
Figure  4.2  shows  the  continuation  result.  The  continuation  result  remains  a  single 
elongate  anomaly  and  any  detection  algorithm  will  likely  treat  it  as  a  single  anomaly. 

In  the  above  case,  the  the  maximum  height  is  2.5  m  and  the  dipole  depth  is  0.5 
m.  The  total  distance  between  the  observation  plane  and  source  is  3.0  m,  which  is  1.5 
times  the  dipole  separation.  In  general,  this  height  depends  on  the  dipole  separation, 
noise  level,  and  acquisition  geometry.  A  higher  signal  to  noise  (SNR)  of  data  would 
increase  this  height.  Given  the  current  quality  of  magnetic  data  in  UXO  application, 
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Figure  4.2.  The  observed  data,  produced  at  an  observation  height  of  3  m  with  two 
dipoles  buried  at  0.5  m  below  the  ground  surface,  is  downward  continued  by  3  m.  The 
result  from  data  at  this  height  do  not  have  any  indication  of  two  separate  anomalies, 
thus  showing  the  limitations  of  the  approach. 
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Figure  4.3.  Comparison  of  the  observations  at  two  different  heights  (2.5  m  and  3.0 
m)  above  the  ground  surface,  respectively,  and  the  corresponding  results  of  stable 
downward  continuation. 
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however,  we  can  infer  that  an  observation  height  comparable  to  the  dipole  separation 
(1.0  to  1.5  times)  is  a  reasonable  limit. 

4.2  Data  Enhancement 

We  turn  our  attention  to  the  application  of  SDC  algorithm  to  field  data  and 
examine  the  potential  enhancement  from  the  continuation  algorithm.  We  apply  the 
method  to  a  small  subset  of  airborne  magnetic  data  acquired  from  helicopter  platform 
over  Kirtland  Air  Force  Base  (Billings  et  al.,  2009).  We  chose  this  subset  because  of 
the  availability  of  ground  magnetic  data  in  the  northern  portion  of  this  area.  The 
overlapping  data  sets  allows  us  to  objectively  assess  the  improvement  in  the  continued 
data. 

Figure  4.4  shows  the  data  set  over  an  area  of  100  m  by  100  m.  The  data  have 
been  gridded  to  an  interval  of  0.25  m  in  both  directions.  The  radial  power  spectrum 
and  its  three-ensemble  representation  with  noise  component  are  displayed  in  Figure 
4.5.  Using  this  power  spectrum  to  obtained  the  corresponding  spectral  weighting  at 
1.0  m  below  the  observation  height,  we  obtain  the  stably  downward  continued  data 
shown  in  Figure  4.6. 

Two  observations  can  be  made.  First,  there  is  significant  striping  in  the  east- 
west  direction.  This  has  been  ascertained  as  the  heading  error  in  the  original  data. 
Although  the  data  had  been  leveled  during  the  post-acquisition  processing,  it  is  clear 
that  the  heading  errors  were  not  completely  removed.  Downward  continuing  the  data 
has  accentuated  this  component  to  the  degree  it  is  visible  now.  The  heading  error 
generally  appears  to  be  relatively  high  frequency  in  the  cross  line  direction  but  long 
wavelengths  in  the  flight  direction.  For  this  reason,  the  stable  downward  continuation 
tends  to  retain  it  in  the  result.  The  identification  of  the  heading  errors  suggests  that 
current  practice  in  leveling  of  airborne  magnetic  data  should  be  improved.  This  is 
especially  important  if  data  are  to  be  used  in  discrimination  process.  Because  of  this 
observation,  our  subsequent  work  with  field  magnetic  data  has  incorporated  a  step 
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Figure  4.4.  A  small  subset  of  Kirtland  airborne  magnetic  data. 

of  micro-leveling  to  treat  the  heading  error.  Performing  downward  continuation  after 
the  required  micro-leveling  on  Kirtland  airborne  data  over  large  areas  no  longer  show 
any  noticeable  heading  errors.  This  will  be  exemplified  in  the  following  chapter. 

The  second  observation  is  that  the  continued  data  are  noticeably  sharper  and 
have  accentuated  higher  frequency  content.  Several  confined  anomalies  are  now  easily 
identified  that  were  not  clearly  visible  in  the  original  data.  The  circles  in  Figure  4.6 
mark  several  such  anomalies.  Each  marked  anomaly  exhibits  the  dipolar  pattern  in 
the  downward  continued  data,  but  they  appeared  in  general  as  simple  highs  in  the 
original  data.  The  dipolar  anomalies  can  be  readily  picked  by  a  reliable  automatic 
anomaly  detection  method. 

We  compare  the  downward  continued  airborne  data  with  the  ground  magnetics 
in  the  same  area  in  Figure  4.7.  This  is  one  of  a  few  areas  that  have  complete  overlap 
between  the  airborne  and  ground  data.  Again  we  have  marked  the  dipolar  anomalies 
revealed  by  the  downward  continuation  and  the  corresponding  anomalies  that  are 
clearly  visible  in  the  ground  data.  Therefore,  it  is  clear  that  the  stable  downward 
continuation  has  performed  as  expected  in  this  case  and  is  able  to  enhance  the  airborne 
magnetic  data. 
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Figure  4.5.  Radial  power  spectrum  of  the  Kirtland  airborne  magnetic  data  shown 
in  Figure  4.4.  The  three  ensembles  correspond,  respectively,  to  the  surface  magnetic 
soil  and  survey  errors,  UXO  at  the  site,  and  a  long-wavelength  component  possibly 
related  to  deeper  geology. 
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Figure  4.6.  Stable  downward  continuation  of  the  Kirtland  airborne  magnetic  data 
shown  in  Figure  4.4.  The  continuation  process  has  reveal  the  presence  of  heading 
errors  that  were  not  completely  removed  by  leveling  of  the  original  data.  This  indi¬ 
cates  the  need  for  further  de-corrugation  of  the  leveled  data.  More  importantly,  the 
continued  data  reveal  clearly  a  number  of  dipolar  anomalies  (marked  by  circles)  that 
were  not  apparent  in  the  original  data  (Figure  4.4). 
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Figure  4.7.  Comparison  between  downward  continued  data  and  ground  observation. 
The  dipolar  anomalies  marked  by  the  circles  in  the  upper  panel  were  revealed  by 
the  downward  continuation  process,  and  they  all  correspond  to  the  same  anomalies 
shown  in  the  ground  magnetic  data  in  the  lower  panel. 

4.3  Summary 

In  this  chapter,  we  have  examined  the  potential  enhancement  of  magnetic  data 
using  synthetic  and  field  data  sets.  The  first  aspect  is  the  enhancement  of  data 
resolution  to  distinguish  between  adjacent  dipolar  anomalies.  We  have  concluded 
that  when  the  observation  height  is  comparable  or  slightly  greater  than  the  lateral 
separation  of  dipolar  sources,  it  is  possible  with  the  current  data  quality  to  separate 
such  anomalies  through  stable  downward  continuation.  This  puts  the  algorithm  well 
within  the  parameter  range  of  current  helicopter-based  airborne  survey  for  UXO 
applications. 

Application  to  a  small  subset  of  data  from  Kirtland  Air  Force  Base  clearly  shows 
that  several  dipolar  anomalies  are  revealed  through  the  downward  continuation  pro¬ 
cess  while  controlling  the  noise  magnification.  Comparison  with  ground  magnetic 
data  show  that  the  enhance  dipolar  anomalies  are  real  and  correspond  to  the  same 
anomalies  present  in  the  ground  data.  Thus,  we  have  demonstrated  the  enhancement 
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of  anomaly  in  actual  data  acquired  for  UXO  applications. 

An  added  benefit  of  the  work  is  the  identification  of  heading  errors  that  exist  in 
the  leveled  airborne  data.  Although  such  subtle  error  is  not  visible  in  the  original 
data,  stable  downward  continuation  is  able  to  identify  it  because  of  its  long- wavelength 
nature  along  flight  lines.  Thus,  further  it  is  necessary  to  carry  out  de-corrugation  of 
such  data  prior  to  further  analysis.  In  the  next  chapter,  we  conclude  the  project  by 
applying  the  developed  package  and  analyzing  its  performance  to  a  series  of  large- 
scale  HcliMag  System  data  collected  as  part  of  the  ESTCP  wide  area  assessment 
(WAA)  demonstration  project  (MM-0741)  at  Kirtland  Air  Force  Base. 
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CHAPTER  5 


FIELD  EXAMPLE 


The  data  analyzed  in  this  chapter  were  acquired  as  a  part  of  ESTCP  wide  area 
assessment  (WAA)  demonstration  project  (MM-0741)  at  the  Kirtland  Air  Force  Base 
west  of  the  city  of  Albuquerque,  New  Mexico.  ESTCP  established  a  6,500  acre  demon¬ 
stration  sub-area  for  the  WAA  Pilot  Program  within  the  15,246  acre  FUDS.  Docu¬ 
mented  ordnance  at  site  includes  M38A2  100-lb  practice  bombs,  M85  100-lb  practice 
bombs,  and  250-lb  general  purpose  HE  bombs. 

The  airborne  magnetic  data  were  acquired  using  the  HcliMag  System  operated 
by  Sky  Research,  Inc.  There  are  13  sensors  spaced  0.75  m  apart  on  the  platform. 
Flight  lines  were  at  7  m  to  provide  a  40%  overlap.  The  altitude  ranged  from  1  to  3 
m.  Data  density  along  track  was  one  reading  per  0.2  m  (Billings  et  ah,  2009). 

We  have  processes  four  different  areas  within  the  survey.  For  brevity,  we  present 
the  result  from  one  area,  which  is  the  Central  South  grid  shown  in  Figure  5.1  and 
covers  approximately  300  m  by  300  m. 

Figure  5.2  displays  the  de-corrugated  data  after  gridding  to  an  interval  of  0.25 
m  in  both  directions.  The  color  bar  corresponds  to  histogram  equalized  color  con¬ 
tour  intervals.  Numerous  localized  anomalies  are  visible.  This  data  set  has  been 
further  processing  using  a  de-corrugation  filter  to  remove  the  residual  heading  error 
as  observed  in  the  preceding  chapter. 

We  have  stably  downward  continued  these  data  by  1  m.  Figure  5.3  shows  the 
associated  data  misfit,  model  objective  function,  and  choice  of  optimal  regularization 
parameter  by  the  L-curve  criterion.  The  optimal  the  algorithm  is  able  to  pick  the 
well-defined  corner  on  the  Tikhonov  curve  (lower  left  panel)  as  indicated  by  the  single 
peak  on  the  curvature.  That  is,  the  noise  estimate  in  this  case  is  quite  reliable. 

The  predicted  data,  which  constitute  the  de-noised  data  at  the  original  observa- 
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Figure  3.  Overview  of  survey  data  collected. 

The  areas  with  opaque  palettes  were  collected  with  the  13  sensor  HeliMag  configuration. 
The  large  semi-opaque  palette  is  the  AMTADS  2005  data  set.  and  the  boxes  outlined  in  red 
represent  areas  that  were  surveyed  with  a  ground-based  towed  magnetometer  array  in  2005. 


Figure  5.1.  Figure  from  MM-0741  Report  and  associated  caption  indicating  the  area 
plan  of  the  WAA  survey  and  the  locations  of  the  HeliMag  surveys.  The  data  set 
presented  in  this  chapter  is  from  the  Central  South  grid. 


tion  surface,  are  shown  in  Figure  5.4.  This  set  of  data  is  a  high  fidelity  representation 
of  the  original  data  and  the  difference  ranges  from  -12  to  12  nT.  At  the  scale  displayed 
in  the  figure,  little  difference  is  visible.  The  downward  continued  data  are  displayed 
in  Figure  5.5,  which  exhibit  much  greater  amplitudes  and  most  anomalies  are  visibly 
sharper.  There  is  no  longer  visible  stripping  in  the  data  due  to  the  heading  errors. 

To  verify  the  data  enhancement  gained  through  the  stable  downward  continua¬ 
tion,  we  next  apply  an  automated  anomaly  detection  algorithm  based  on  the  extended 
Euler  deconvolution  (Davis  and  Li,  2010).  The  algorithm  was  developed  initially  with 
funding  support  from  ERDC  and  SERDP  projects  MM-1414  and  MM-1638.  The  au¬ 
tomated  detection  solves  for  the  structural  index  (SI)  of  each  anomaly  and  treats  those 
having  an  SI  greater  than  a  threshold  as  potential  UXO  (e.g.  near  3)  for  metallic 
targets.  A  statistical  analysis  based  on  the  amplitude-approximated  source  strength 
discard  spurious  solutions  and  the  remainders  are  kept  as  target  picks.  A  cluster  algo¬ 
rithm  based  on  proximity  of  these  picks  are  then  applied  to  group  duplicate  solutions 
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Figure  5.2.  Airborne  magnetic  data  acquired  by  HeliMag  system  of  Sky  Research, 
Inc.  The  data  are  located  over  the  Central  South  grid  as  indicated  in  Figure  5.1 
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Figure  5.3.  Data  misfit,  model  objective  function,  and  associated  L-curve  picking 
of  the  optimal  regularization  parameter  for  the  stable  downward  continuation  of  the 
Kirtland  data  in  Figure  5.2 


into  final  picks.  For  the  current  data  set,  the  SI  threshold  was  chosen  as  2.5  and  the 
clustering  distance  was  1.0  m. 

Figure  5.6(a)  is  the  result  of  automated  detection  from  the  original  de-corrugated 
airborne  data,  and  Figure  5.6(b)  shows  the  result  from  the  downward  continued  data. 
Both  plots  have  their  respective  data  set  as  color  images  in  the  background  and  white 
circles  denoting  the  target  picks.  The  images  cover  the  same  area  as  the  preceding 
data  figures. 

At  first  look,  the  two  sets  of  results  appear  similar.  However,  closer  inspection 
shows  that  there  are  two  major  differences.  First,  the  results  from  SDC  data  have 
more  individual  anomalies  picked,  indicating  the  continuation  indeed  accentuated 
anomalies  so  that  the  detection  algorithm  is  able  to  identify  them.  This  is  the  same 
enhancement  as  illustrated  through  the  small  data  set  in  the  preceding  chapter,  where 
several  broad  anomalies  are  sharpened  up  to  show  their  correspondence  with  the 
ground  data.  Secondly,  the  majority  of  the  anomalies  detected  in  the  SDC  data  have 
single  circles  marking  each  of  them,  whereas  the  same  anomalies  in  the  original  data 


36 


Predicted  Data 


3894000 


3893900 


3893800 


145.89 
6.20 
2.44 
1.27 
0.67 
0.29 
0.01 
-0.22 
-0.42 
-0.60 
-0.76 
-0.92 
-1.09 
-1.26 
-1.45 
-1.65 
-1.89 
-2  22 
-2.70 
-3.61 
-34.76 


333200 


333300 


333400 


Easting  (m) 


Figure  5.4.  Predicted  data  produced  during  the  stable  downward  continuation  of  of 
the  data  shown  in  Figure  5.2.  This  would  be  the  de-noised  data.  They  maintain  the 
fidelity  to  the  original  data,  and  there  is  little  visual  indication  of  large  differences. 
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Figure  5.5.  Stably  downward  continued  data  of  the  Central  South  grid  shown  in 
Figure  5.2 
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tend  to  have  multiple  circles  that  are  slightly  shifted  from  each  other.  That  is,  each 
anomaly  in  the  original  data  has  produced  multiple  picks.  This  occurs  when  the 
anomalies  are  broad  and  affected  by  noise.  The  enhanced  anomalies  through  SDC 
are  sharper  due  to  the  increased  resolution  and  are  less  affected  by  noise.  Therefore, 
the  detection  algorithm  is  able  to  produce  much  cleaner  picks. 

To  illustrate  these  observations,  we  zoom  in  on  two  smaller  areas  within  this 
grid  to  show  more  details.  Figure  5.7  shows  the  comparison  between  detection  results 
from  the  original  data  (a)  and  SDC  data  (b)  in  the  northeast  corner  of  the  grid.  This 
area  has  sparsely  distributed  dipolar  anomalies.  Primary  improvement  is  indicated 
by  mostly  single  picks  for  each  anomaly  in  SDC.  There  are  several  occurrences  where 
increased  data  resolution  resulted  in  more  picks,  for  example,  towards  the  upper  left 
corner,  center  of  the  map,  and  lower  left  corner. 

Figure  5.8  shows  a  similar  comparison  between  detection  results  from  the  original 
data  (a)  and  SDC  data  (b)  in  the  southwest  corner  of  the  grid.  Many  more  dipolar 
anomalies  are  closely  located,  especially  towards  the  west  edge  of  the  map.  Because  of 
the  interference  between  adjacent  anomalies  in  the  airborne  data,  many  of  them  are 
missed  by  the  automated  Euler  detection  algorithm.  However,  the  stable  downward 
continuation  has  significantly  increased  the  resolution  and  enabled  the  separation  of 
these  anomalies.  Consequently,  the  Euler  detection  is  able  to  capture  nearly  all  of 
the  dipolar  anomalies,  ft  only  missed  three  large-scale  anomalies,  each  spanning  over 
10  to  15  m,  which  is  much  wider  than  the  maximum  window  size  of  5.25  m  used 
in  the  detection.  The  difference  displayed  in  this  figure  clearly  demonstrates  the 
enhancement  which  the  stable  downward  continuation  algorithm  is  able  to  produce. 
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Figure  5.6.  (a)  Automated  Euler  detection  results  from  the  original  data  in  Figure 
5.2.  (b)  Automated  detection  from  SDC  data  shown  in  Figure  5.5  by  using  exactly 
the  same  parameters. 
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Figure  5.7.  Detailed  plots  of  northeast  corner  of  the  gird,  (a)  Automated  Euler 
detection  results  from  the  original  data,  (b)  Automated  detection  from  SDC  data. 
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Figure  5.8.  Detailed  plots  of  southwest  corner  of  the  gird,  (a)  Automated  Euler 
detection  results  from  the  original  data,  (b)  Automated  detection  from  SDC  data. 
Euler  detection  is  able  to  capture  nearly  all  of  the  dipolar  anomalies  from  the  SDC 
data. 
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CHAPTER  6 


SUMMARY 


Successful  UXO  discrimination  is  essential  for  minimizing  the  costs  of  ordnance 
clearance.  Discrimination  algorithms  depend  on  accurate  and  reliable  sensors  for  dif¬ 
ferent  types  of  data  and  the  availability  of  high  quality  data  these  sensors  provide. 
Much  of  the  attention  in  UXO  geophysics  research  has  been  devoted  to  the  devel¬ 
opment  of  new  sensors  and  new  algorithms  for  data  inversion  and  discrimination. 
However,  a  third  and  equally  important  aspect  is  the  pre-processing  of  data  and  the 
extraction  of  data  with  high  signal-to-noise  ratio  (SNR).  Such  algorithms  provide  a 
complementary  component  to  newly-developed  sensors  and  strengthen  the  foundation 
for  advanced  discrimination  criteria. 

The  magnetic  method  is  one  of  the  most  effective  tools  for  UXO  detection  and 
discrimination.  A  great  amount  of  data  is  being  acquired  using  a  variety  of  man- 
portable,  vehicle-towed,  and  airborne  platforms.  The  data  are  used  in  both  discrim¬ 
ination  for  active  clearance  and  wide  area  assessment  (WAA)  (e.g.,  ESTOP  Project 
UX-0533).  In  particular,  recent  developments  have  enabled  the  recovery  of  dipole 
or  ellipsoidal  parameters  through  inversion  of  magnetic  data  (Billings,  2004).  The 
reliability  of  these  algorithms  depends  on  well-defined  anomalies  associated  with  in¬ 
dividual  targets  and  a  quantitative  characterization  of  data  noise. 

The  signal  content  and  relative  noise  characteristics  in  a  magnetic  data  set  are 
strongly  affected  by  the  observation  height.  It  is  well  known  that  the  magnetic 
anomaly  due  to  a  dipolar  source  decays  with  distance  cubed  and  the  anomaly  spreads 
out  spatially  in  proportion  to  the  distance  as  well.  According  to  (Telford  et  al.,  1990), 
for  example,  the  width  of  a  magnetic  dipole  anomaly  at  one  half  of  its  maximum  is 
approximately  equal  to  the  source  depth.  This  means  that  the  data  resolution  defined 
as  the  distance  between  separable  anomalies  decrease  linearly  with  increasing  obser- 
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vation  height  and  the  relative  noise  in  the  data  increase  with  the  observation  height. 
The  ideal  data  are  those  acquired  with  a  sensor  at  zero  height  above  the  ground. 
However,  because  of  the  intrinsic  limitations  imposed  on  field  data  acquisition,  all 
previously-mentioned  platforms  collect  data  at  some  height  above  the  ground  sur¬ 
face.  For  example,  the  magnetic  sensor  height  in  a  typical  man-portable  system  or  a 
vehicle-towed  system  (e.g.,  MTADS)  may  be  between  0.15  to  0.5  m  above  the  ground, 
and  a  helicopter-based  survey  would  be  flown  at  a  height  about  1.5  to  2.0  m  above 
the  ground  (Doll  et  al.,  2006).  Thus,  there  is  a  need  for  enhancing  magnetic  data  to 
alleviate  the  effects  of  observation  height  and  to  estimate  the  noise  characteristics. 
This  SEED  project  was  established  to  address  this  specific  aspect  of  magnetic  data 
by  developing  and  applying  a  stable  downward  continuation  method. 

6.1  Proposed  Research  Accomplishments 

During  the  first  year  of  the  SEED  project  the  main  tasks  were  completed  with 
a  working  SDC  algorithm  developed  and  tested  on  synthetic  and  field  data.  The 
algorithm,  presented  in  Chapter  2,  formulates  downward  continuation  as  an  inverse 
problem  using  Tikhonov  regularization  and  has  the  flexibility  of  incorporating  the 
expected  power  spectrum  of  UXO  anomalies.  The  degree  of  regularization  for  this 
formulation  was  estimated  automatically  using  the  well-established  methods  in  lin¬ 
ear  inverse  problems.  Applications  during  the  first  year  showed  that  the  algorithm 
reliably  estimates  the  noise  in  UXO  data  and  reconstructs  the  magnetic  anomaly  at 
ground  surface  within  the  limitation  imposed  by  the  noise.  The  reconstructed  field  at 
the  ground  surface  exhibits  significant  enhancement  compared  to  the  original  data. 

6.2  Expanded  Research  Accomplishment 

The  second  year  of  the  project  significantly  advanced  upon  the  original  proposed 
research  by  incorporating  an  initial  set  of  inversion  algorithms  to  estimate  requisite 
ensemble  depths  within  the  data’s  radial  power  spectrum  for  improved  SDC.  These 
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developments  are  presented  in  Chapter  3  of  the  report.  In  particular,  development 
of  recursive  Quenched  Simulated  Annealing  (QSA)  from  SERDP  Project  MM-1638 
has  been  adapted  for  this  purpose.  The  QSA  algorithm  has  demonstrated  a  reliable 
tool  for  modeling  the  redial  power  spectrum  of  UXO  magnetic  data  using  non-linear 
parametric  inversion.  As  a  result,  we  can  now  automatically  define  a  model  objective 
function  for  SDC  to  quantify  the  conformity  of  the  enhanced  data  with  expected 
spectral  properties.  In  addition,  we  have  been  able  to  validate  the  SDC  enhancement 
package  over  the  last  year  by  evaluating  its  performance  with  the  latest  detection 
technology  based  on  the  extended  Euler  deconvolution. 


6.3  Future  Research  Recommendations 

This  project  has  demonstrated  the  validity  of  enhancing  magnetic  data  using 
SDC  and  it  is  currently  a  well  developed  package  for  application  in  active  UXO  re¬ 
mediation  efforts.  However,  there  are  still  unanswered  research  questions  necessary 
to  complete  development  and  understanding  of  the  approach.  The  first  question  is 
understanding  the  effects  and  limitations  of  unavoidable  variation  in  sensor  height 
due  to  platform  motion.  This  is  particularly  relevant  for  airborne  systems  imple¬ 
mented  during  WAA  deployment.  The  second  area  is  building  an  understanding  of 
scale  for  practical  implementation  of  SDC.  We  have  identified  through  this  project 
that  magnetic  data  at  UXO  sites  will  require  sectioning  into  sub-grids  for  proper 
data  enhancement.  A  method  must  be  developed  that  allows  us  to  objectively  iden¬ 
tify  sub-regions  within  the  remediation  site  such  that  the  data  demonstrate  similar 
spectral  and  statistical  properties  for  effective  SDC  application.  Finally,  SDC  method 
has  demonstrated  the  potential,  through  ensemble  based  reconstruction  of  a  model 
objective  function,  the  ability  to  selectively  downward  continue  data  for  enhanced 
data  analysis.  The  natural  result  is  that  it  is  possible  to  expand  the  method  into 
the  realm  of  ’geology/background  separation’  algorithms  needed  in  UXO  remediation 
efforts,  particularly  for  difficult  data  sets  and  environments.  For  example,  it  is  likely 
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that  SDC  can  contribute  to  the  enhancement  of  difficult  magnetic  data  at  sites  such 
as  San  Lnis  Obispo  in  California. 
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